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Abstract. Infiltration of anomalously diffusing particles from one material to 
another through a biased interface is studied using continuous time random walk and 
Levy walk approaches. Subdiffusion in both systems may lead to a net drift from one 
material to another (e.g. (x(t)) > 0) even if particles eventually flow in the opposite 
direction (e.g. number of particles in x > approaches zero). A weaker paradox is 
found for a symmetric interface: a flow of particles is observed while the net drift is 
zero. For a subdiffusive sample coupled to a superdiffusive system we calculate the 
average occupation fractions and the scaling of the particles distribution. We find a 
net drift in this system, which is always directed to the superdiffusive material, while 
the particles flow to the material with smaller sub or superdiffusion exponent. We 
report the exponents of the first passage times distribution of Levy walks, which are 
needed for the calculation of anomalous infiltration. 
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1. Introduction 

Infiltration of diffusing particles from material A to material B through some interface is 
a widely investigated process. In recent years much focus was diverted to the problems 
when the diffusion in one material or in both is anomalous, namely (x 2 (t)) oc t a and 
a 7^ 1 [U El [3] . This behavior is important in numerous applications such as infiltration 
of water into porous soil [HE], contaminant diffusion [6J, [T| , moisture ingress in zeolites 
[8] or in fired clay ceramics [9j \TU\ , diffusion of sugar through a membrane in a gel solvent 
[TT] . and polymer translocation through a nanopore [T2J [T3], [H] (see also [T5J [TB] [T7]). 

We consider two semi-infinite materials located in x < and x > 0, where particles 
exhibit anomalous sub or superdiflusion. We find several peculiar behaviors unique to 
anomalous infiltration: in the case of a composition of two subdiffusive systems the 
infiltration of particles from material x < to x > leads to a net drift (x(t)). Such 
a drift increases slower than t 1 ^ 2 as expected from unbiased diffusion processes, still it 
may yield a net sub-current j — d (x) /dt, which vanishes as t — > oo. We also show that 
in some cases the flow of particles is opposite to the drift. In fact we find a situation 
when asymptotically all the particles are say in sample x > but the average drift (x(t)) 
is oppositely directed (x(t)) < 0. This seemingly paradoxical behavior is explained in 
the text. Secondly, if materials x < and x > have different diffusive properties, 
in the long time limit all particles will be accumulated in the material with slower 
diffusion, which will act as a trap producing a flow of particles from one material to 
another. Interestingly, in the long time limit the drift depends only on the properties 
of the slower medium. This is a surprising result since (x(t)) can be very far from the 
interface, deep in the faster sample, but still be independent of the properties of that 
region. 

A second model we study is a subdiffusive material (for example in x < 0) coupled 
to superdiffusive sample (in x > 0). For super diffusion motion we consider a Levy walk 
model [HI [TU [20] . To analyze this infiltration problem we need the distribution of the 
first passage times (FPT) [21J for a Levy walk on a semi-axes, which are reported here 
for the first time (see J22j [23[ [231 [2H] for other works on anomalous first passage time 
problem). Using the FPT density in x > and x < 0, we calculate the average of 
occupation fractions and find the scaling of the particles distribution. For a subdiffusive 
system coupled to a superdiffusive material a net drift is found even for unbiased motion 
on the boundary. This drift is always directed to the superdiffusive material, while the 
particles flow to the material with longer sticking times. 

Although phenomena such as drift against the flow and flow without the drift 
are known for systems with normal diffusion where they are generated by geometrical 
constraints, or by thermal or external field inhomogeneities [26l [27[ |28[ [29] . in our case 
these effects are only due to the anomalous nature of diffusion and are not present for 
normal diffusion. These phenomena are explained by the competition of the diffusion 
processes which are slower or faster than normal spreading. Part of the results were 
shortly summarized in [30J. 
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2. Two Coupled Subdiffusive Systems 

We consider two coupled subdiffusive materials using the continuous time random walk 
(CTRW) model as the underlying process [31]. In the continuum limit the fractional 
diffusion equations in materials x < and x > describe the dynamics [2} [32], [33], [34] [35] 

d - E ^T 1 = oDt a+ K^P(x,t), x>0, (1) 
where the Riemann-Liouville operator oD t ~ a is defined as [2] [36] 

oA p{Xit) -f(a)diJo dt (t-t>y- a} ( ) 

and 0<o~<l,0<o + <l. Constants are anomalous diffusion coefficients 

with units [m 2 /sec a ] and [m 2 /sec a+ ], respectively. As well known, the fractional 
diffusion equation (JTJ) with a - = a + = a and K~ = K + = K yields for particles 
starting on the origin (x 2 {t)) = 2Kt a /T{l + a) [321 [331 (Ml 135] • To solve equation CQ) 
we determine the boundary conditions for this equation starting with a random walk 
picture. 



2.1. CTRW model: Drift 

We consider a CTRW on a one dimensional lattice (see figured]) with the lattice spacing 
a, which in the continuum limit will be made small. For lattice points i<0a particle 
has the probability 1/2 to jump to one of its nearest neighbors. Waiting times on each 
lattice point are independent identically distributed random variables with a common 
PDF if)~(r). For x > a similar unbiased random walk takes place with a waiting time 
PDF ^ + (r). On the lattice point x = (the boundary) a particle has the probability to 
jump right q + or left q~ = 1 — q + and the waiting times are exponentially distributed 
with a rate Rq. Thus, a particle starting on the origin will jump say to the right (with 
probability q + ) after waiting an average time 1/Rq, then on the lattice point x = a, 
it will wait for time r drawn from iP + (t), and then with probability 1/2 will jump to 
the left or right. The waiting times have power law distributions ip~(r) oc r _1 ~ a and 
oc r _1_a+ , as r — > oo. More specifically, using Tauberian theorem the Laplace 
transform r — > s of the waiting time PDF behaves like 

4>-(s) oc 1 - B- S a ~, ip + (s) oc 1 - B + s a+ , (3) 

when s — > corresponding to r — > oo. Throughout the paper we denote the Laplace 
transform as f(s) = / °° dt e~ st f(t). The anomalous diffusion coefficients are given by 
K~ = lim a 2^ 0i £-^ a 2 /2B~ and K + = \im a 2^ B +^ a 2 /2B + [35]. Our results are not 
changed if on x = the waiting times are power law distributed like ip~ or V' + instead 
of exponential. 
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Figure 1. Illustration of the random walk model. The times which particles spend at 
sites in the region x < and x > are distributed by the waiting time PDFs and 
respectively. The waiting times at x = are exponentially distributed. 



Using the CTRW model we find the drift in the long time limit in the following 
way: From the start of the process at t — until time t the particle made N jumps, 
where N is a random variable. Hence the position of the particle at time t is 

N 

x = J2^i, (4) 

i=l 

if at t = the particle is on the origin. From the model assumptions Sxi is equal +a 
or —a and it describes the zth jump in the sequence. The jump length 5xi satisfies 
(Sxi) = if the particle is not on the origin since then the probability of jumping left 
or right is equal 1/2. Therefore 

(x(t))=a(q+-q-)(n g (t)), (5) 

where (n z (t)) is the average number of times the particle visited the origin, which 
is calculated in the Appendix A using the CTRW model. In particular (n z (t)) is 
determined by the first passage time PDFs in samples x > and x < since these 
first passage times determine the number of visits to the origin (see details in Appendix 
A). To derive equation (JSJ) we have used the fact that on the origin the average step size 
is a(q + — q~). In [37] we related q + , q~ to energy gap between material A and B using 
detailed balance condition. When a~ = a + = a 



I u\\ V + ~ % VK K+ /2 

(x(t)) oc — ; 7= t 1 . 6 

V w/ r(l + a/2) q-y/W + q+VlT KJ 

For the case a~ < a + , we get 

( x(t )) oc q+ - q - ^ e~' 2 (?) 

{X[t))CC q- r(l + a-/2) ' [) 
which agrees well with simulations in figure[2j The sign of the drift, i.e. its directionality, 

is determined by the sign of q + — q~ , and (x(t)) = if q + = q~. Equation (JTJ) shows 
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Figure 2. The drift (x(t)) (open circles) and (x 2 ) (filled circles) calculated numerically 
for the CTRW model with a+ = 0.75, K+ = 0.138, or = 0.3, K~ = 0.385 and 
q + = 0.7. Dashed and dashed-dotted lines represent long time asymptotic behavior 
described by equations ([7]) and (|26l) . respectively. 

that in the long time limit the drift depends only on one diffusion constant in sample 
x < (i.e. K~) and grows in time with the exponent of the slower medium (i.e. or). 
This is a surprising result since (x(t)) can be very far from the interface, deep in the 
faster sample x > 0, but still be independent of the properties of that region a + , K + . 
The exponent of the drift (x(t)) cx t a I 2 is the exponent of the slower medium, which is 
clearly related to the power law distribution of first passage times in the slower medium 
(f>~(t) oc / 2 ) [22J (see Appendix A). One can show that equation (JTj) is valid for 

times 

1 « J^ q - 
V K~ q + 



2.2. CTRW model: Statistics of occupation times 

Now we consider the distribution of occupation times in the material x < or x > 0. 
Let ft{t~) be the PDF of the total time t~ a walker stays in the material x < and t 
is the measurement time. Similarly, t + is the total time a walker stays in the material 
x > 0. The double Laplace transform of /*(*"), f s {u) = / °° dte~ st / °° dt~ e~ ut ~ f t (t~), 
reads (the derivation is given in Appendix B) 

q-(s + Mj^/v^ + g+s^/v^ 
For a~ = a + equation (JS} reduces to the Lamperti PDF [38], [39] (see Appendix B). 
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Let us assume a < a + . Expanding equation (jSJ) in u, we get 

r(s))(x ^ , (9) 

where TZ(s) is given by 



a + \/K- s a /2 

R "» = q rim^ (10) 

In the long time limit t — > oo (s — > 0) equation fl9]) gives (after inverting the Laplace 
transform) 



q + {K 1 t 1 



a 1 —a 



t~{t)) (xt-—\ —— ^— . (11) 

w/ g-v^ + r(2-^^) 

Now it is straightforward to get the average occupation fraction in sample x < 




OL — a 1 



= ^ - 1 - g-A/^r — ■ (12) 

t q-y K+ T(2 + ^- p,+ 1 

Since a~ < a + , the second term in this equation vanishes as t — > oo and V~(t) — > 1, 
which indicates that in the long time limit all particles will be located in the region 
x < (see figure [3]), a result valid for any q~,q + < 1. As might be expected, the slow 
domain ct~ < a + serves as a perfect trap: all particles flow to the slower domain. The 
usual normalization condition 

V~(t)+V + (t) = l, (13) 

gives V + (t) in sample x > 0. 

2. 3. Boundary conditions and solution of equation (Q]) 

We now wish to find P(x, t) i.e. solve equation ([1]). To obtain the solution of a standard 
or fractional diffusion equation, one has to know the mathematical boundary conditions 
between sample x < and sample x > 0. One boundary condition is well known 
and needs no further discussion: the probability current must be balanced so that 
normalization is conserved (conservation of number of particles), see some details below. 
Previous works assumed in addition the second boundary condition which demand the 
constant ratio of the concentrations at two sides of the boundary located at x = 0, 
namely P(x,t)\ x=0 - = kP(x, t)\ x=0 +, where k was assumed to be equal to k = 1 [30] or 
k = const. [UJ. For normal diffusion in samples x < and x > n was derived using a 
normal random walk theory [12]. A generalization of the boundary conditions for sub- 
diffusion with unbiased boundary, was considered in [T71 14_3] (see also [35j H3J H3, 06] ) . 
Solution of equation (pQ) in Laplace space reads 

s^/^exp 

P(x, s) = C + (s) J -9(x) + 



,0-/2-1™ _\x\s 



s ' exp 



+C'(s) n ^ ' [1 " 0{x)\ , (14) 
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Figure 3. The occupation fraction in x < 0, V~(t) = f_ dx P(x,t), calculated 
numerically for the CTRW model (open circles). Parameters are the same as in figure 
[21 The dashed-dotted line given by equation (TT2j) describes how V~ approach its limit 
V~ — s- 1 (dashed line). The dotted line corresponds to V~(t) calculated using equation 
(|C.1|) . Notice that all the particles flow to the left, however (x(t)) > namely particles 
drift to the right (see figure [2]). 



where C~(s) and C + (s) are functions soon to be determined. Here P(x,0) = 5(x) is 
used for the initial condition. The conservation of probability: f_™ dxP(x, s) = 1/s, 
gives C + (s) + C~(s) = 2. Using equation fTH]) we find 

_ s - a++lK+ dP(x,s) 



-i K - dP(x,s) 
dx 



\x=0~ 



dx 



1 2=0-1 



(15) 



For a < 1 and a + < 1, Equation (fTBI) represents the continuity of fractional probability 
flow at the boundary 



J+( X = o+,t) - J-(x = 0~,t) = 8(t). 



(16) 



The fractional probability flow in this case is the generalization of the usual definition 
[34J, for example for x < 



J-(x,t) = -K- Dl 



dx 



(17) 



Therefore the fractional equation equation ([T]) can be written as 

dP(x,t)/dt + dJ-/dx = 0, (18) 

for x < and similarly for x > 0. Thus, equation (1151) is nothing else but the Laplace 
transform of the condition of continuity of the fractional probability current at the 
boundary J~(x — 0~,t) = J + (x = + ,t). 
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Figure 4. PDF of particle's position P(x, t) calculated numerically by CTRW model 
(dotted lines) perfectly agrees with analytical theory (dashed lines) found inverting 
equations (fT4"]) and (|20|) to the time domain. The figure illustrates that the majority 
of the particles are found in the slow domain x < 0, however the tail of the PDF 
extends deeply into the fast domain. Thus even if eventually all the particles will be 
accumulated in the slow domain on the left, (x) will be located in the fast domain on 
the right. We used a~ = 0.3, K~ = 0.385, a+ = 0.8, K+ = 0.108 and q+ = 0.7 at 
time t = 10 6 (solid line) and t = 10 8 (dotted line). For t = 10 8 almost all particles are 
in x < 0, V~ ~ 0.93, but the drift is positive (x) = 14.2. 



To derive the second boundary condition we calculate the first moment, (x(s)) = 
fZ^ dx x P(x, s) (see Appendix C for the alternative derivation of a second boundary 
condition). Using equation ( fl4l) 



(x( s )) = — ( Vk+C + (s)s-^ - ^lFc-{s)s~^ ) . (19) 



We require that equation (fl9l) be equal to the Laplace transform of equation (jij), which 
was calculated from the CTRW model. For a~ < a + , equations ffl9l) and (j7]) yields 
when s — > 




C + {s)~^— a-*- , C-(s) = 2-C+(s). (20) 

Using equations (1201) and jT]), we derive the second boundary condition 

q +K~s- a 'P(x,s)\ x=0 - = q-K + s~ a+ P(x,s)\ x=0+ . (21) 

Equation (1211) shows that generally the PDF at the boundary is not continuous, similar 
to the normal diffusion case [12]. Such a jump in the PDF on the origin is shown in 
figure HI Equation (1211) also shows that the scaling of the solution is more complex than 
in time- fractional diffusion equation with one diffusion exponent. Using equations ( fL~ij) 
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Figure 5. Same as in figure U presented in scaled form given by equations 



and f )20|) we find the scaling of the solution for a < a + 

{ t a ~ /2 P(x,t), forx<0 
G a - a+ {0 = I . " (22) 
' \t a+ - a -' 2 P(x,t), forx>0, V } 

and 

f \x\/t a ~ /2 , for x < 
£=< +l (23) 
\x/t a+ '\ forx>0. 

Numerically calculated scaled PDF is shown in figure |5] and exhibits data collapse. 

It is straightforward to calculate moments using the analytical expression for the 
propagator equations ( 1T4"1) and ( [201 . For the mean we obtain 



/ / ^ (q + - q-)VK-K+ 
(x(s)) = == (24) 



Expanding equation f[24l in s we get the result, which coincides with one calculated 
from the CTRW model (see equation (jBJ)). For the second moment we get the following 
expression 

2VK-K+ (q-y/K = 8- a ~f 2 + q+VK+s~ a+ / 2 



(* 2 (s)) = 7 _,„ . ^= L , (25) 



which in the long time limit t — y oo (s — y 0) gives for a - < a 



q + 2VK-K+6 a ' +a+ y 2 



x {t) ) ~ F r ( i + ( a- + «+)/2) - (26) 

This result is in agreement with numerical simulations of the CTRW model (see figure 
W). 
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2-4- Drift against the flow and flow without the drift 

Consider the situation where a~ < a + , so we call the domain x < the "slow" medium. 
From equation ( fT2i) . the probability to be in the slow medium in the long time limit is 

p-(t)-tl, t^oo. (27) 

We also calculated the mean drift in this case (see equation (JTj)). Thus, as mentioned, 
independently of the details of the model all particles flow into the slower medium, 
which absorbs them in the long time limit. However, at the same time if q + > q~ , 
the drift (x) > is positive and increasing with time. Namely, (x) is located in the 
"fast" medium even though all particles eventually accumulate in the slow medium. As 
mentioned, while the dynamics in the faster domain x > is clearly important (since 
(x) may be in that domain) the mean (x) does not depend on the diffusion constant 
K + of that medium, neither on the anomalous diffusion exponent a + . 

An explanation of this paradox is as follows: Although the region with smaller a 
will accumulate more and more particles in the long time limit, at finite time there will 
be always some particles in the opposite region where a is larger. As shown in figure IU 
these particles are moving more freely and travel far away from the interface which will 
compensate the accumulation of particles in the region with smaller a. In other words 
while V + = 1 — V~ = P(x,t)dx — > 0, which naively implies lim^oo P(x,t) = for 
x > (since P(x,t) > 0), still J °° x P(x,t)dx does not approach zero. 

3. Coupled Sub and Superdiffusive systems 

Now we consider a composition of subdiffusive material in one region (for example x < 0) 
and a material with super diffusion in the other region (x > 0). Subdiffusion is modeled 
by the CTRW on a lattice. As before, a particle has the probability 1/2 to jump 
to one of its nearest neighbors. Waiting times on each lattice point are independent 
identically distributed random variables with a common PDF ip(r) oc r _1_a as r — > oc 
with < a < 1. Thus, in our composite system a particle starting on the origin will 
jump say to the right after waiting a random time drawn from the distribution ip(r), 
and performs a super- diffusive walk in x > until it returns and crosses the boundary 
x = 0. Then, a particle performs CTRW in x < until it returns the boundary and so 
on. At the boundary we consider equal probabilities to go left or right q + = q~ = 1/2. 

For superdiffusion we consider the Levy walk model, which corresponds to the 
spatiotemporally coupled version of the CTRW [HI [20]. The waiting time and jump 
length PDFs are no longer decoupled but appear as ip(x,t) = \(x)p(t\x). We consider 
the coupling in the form p(t\x) = |5(|x| — vt), where v is a velocity and the PDF of 
jump length X(x) oc <t 7 /|x| 1+7 with < 7 < 2. In what follows we consider v = 1. Since 
the velocity v is finite it penalizes long jumps such that the overall process attains a 
finite variance (in contrast to infinite variance of Levy flights (x 2 (t)) = 00) [T9| [20] . For 
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a Levy walk (i.e. without coupling to a subdiffusive system) 
, \ f t 3 ~ J , for 1 < 7 < 2, 

x *)* +2 f n ~ (28) 

7 If, for < 7 < 1. 

For 7 = 2 the Levy walk converges to Gaussian process with (x 2 (t)) oc t. 

Our aim is to calculate the occupation fractions, the PDF and the drift for the 
introduced coupled sub and super diffusive systems. For this we need to know the first 
passage time (FPT) density of the Levy walks. Let us consider the FPT distribution 
for the Levy walks on a semi axis. Since for 7 — > 2 the Levy walk converges to Gaussian 
process, the PDF of the FPTs for the Levy walk with 7 = 2 should be 0(r) oc r~ 3 / 2 |21j. 
For 1 < 7 < 2 we find that the PDF of the FPTs for the Levy walk to be independent 
of 7 while for < 7 < 1 the PDF of the FPTs depends on 7 

f r~ 3/2 , for 1 < 7 < 2, 

^ a -1-7/2 f n ;/ ( 29 ) 

[ r 11 , for < 7 < 1. 

These results can be deduced from the exponents the FPT densities of the subdiffusive 
CTRW. Coupling a CTRW system with a Levy walk system we expect that occupation 
fractions in both systems will attain finite values only if 7 = a. For that to happen 
FPT in both systems must behave similarly. Hence for a Levy walk with < 7 < 1 
the first passage time PDF is <j>{r) oc r -1-7 / 2 in the super diffusive system corresponding 
to <p(r) oc t -1- "/ 2 in the subdiffusive system. For 1 < 7 < 2 we expect <j>{r) oc r~ 3//2 
since in that regime the coupled system exhibits normal behavior (1 < 7 < 2 implies 
normal first passage times). Numerically calculated FPT densities for the Levy walks 
are shown in figure [6] and they are compatible with equation ( l29l) . 



3. 1 . Occupation fractions 

Using the FPT density we can now calculate the distribution of occupation times in 
the coupled sub and super diffusive systems. For this we use the method described in 
section 2.2 (see also Appendix B). The PDF of occupation times is determined by the 
first passage time PDFs in x > and x < 0. For the CTRW in the material x < the 
first passage PDF is given by 4>~(t) oc r~( 1+a / 2 ) with < a < 1. For Levy walk in the 
material x > the first passage PDF is as we just have shown + (r) oc r~ 3 / 2 if 1 < 7 < 2 
and 4> + {r) oc r~( 1+7 / 2 ) if < 7 < 1. So, depending on the values of a and 7 three cases 
can be classified: Case (I) 1 < 7 < 2 and < a < 1, so a < 7. In this case the average 
occupation fraction in the material x < behaves as V~(t) oc 1 — t _ ( 1_Q )/ 2 — > 1 as 
t — > 00 and in the long time limit almost all particles will be in the material x < 0. 
Accordingly, since V + {t) + V~{t) = 1 

v + (t) oc r {l - a)/2 ->■ 0, t 00. (30) 

For < 7 < 1 there are two cases: case (II) < a < 7 < 1 and case (III) 
< 7 < a < 1. When a < 7 the average occupation fractions in x < and x > 
behave as V~(t) oc 1 - r^V 2 -> 1 and V + (t) oc r^-")/ 2 -»■ as t ->■ 00. For a > 7, 
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Figure 6. PDF of first passage time 4>(t) calculated numerically for the Levy walk 
on a semi axis for (a) 7 = 0.2, (b) 7 = 0.5, (c) 7 = 1.7, (d) 7 = 2. Blue dashed 
lines are proportional to oc r -1-7 / 2 and red dotted lines are proportional to oc t -3 / 2 . 
Simulations agree with our theoretical prediction equation (|29l) . 



V~(t) oc r( a -r)/ 2 -> and correspondingly V + {t) oc 1 - H a - 7)/2 -> 1 as t -> 00. 
Collecting results we write 

f 1, (I): 1 < 7 < 2, < a < 1, 
P"(t) - 0, (II): < 7 < a < 1, (31) 
[ 1, (III): < a < 7 < 1, 

and V~(t) +V + (t) = 1. This result is very natural, wherever we find the largest sticking 
time the occupation fraction in that system will be 1. Only when a = 7, V~ and V + 
are not trivial in the long time limit. 



3.2. Scaling of the PDF 

Using the average occupation fractions we can find the scaling of the particles density. 
First we consider case (I), namely 1 < 7 < 2 and < a < 1. For the coupled sub and 
super diffusive system we are looking for the PDF in x < in the form (see equation 

flU) 



^ W 2 
"exp I— 1 \„ , 

P(x,s) oc y-(s) ^ Vi ^ Q 7 , x<0, (32) 



7 Kn 
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Figure 7. Scaled PDF G Qj7 (£) for the coupled sub and superdiffusive system defined 
in (|39j) calculated numerically at t = fO 3 and t = fO 4 (solid and dashed lines) with 
a = 0.7 and 7 = 1.2. The scaling variable £ is defined in equation (|40|) . The peaks are 
a result of ballistic paths, i.e. path with no turn overs. 



where K a is the fractional subdiffusion coefficient. Integrating equation ( )32|) we find the 
average occupation fraction in x < 0, V~{s) = J oo dxP(x, s) oc Y~(s)/s. However, as 
we have just shown, for case (I) the average occupation fraction V~(t) — > 1 as t — > 00 
or V~(s) oc 1/s as s — > 0, which yields Y~(s) oc 1. For x > 0, similar to equation ( 132|) . 
we are looking for PDF, which in Fourier-Laplace space has the form 

P(k,s) = Y + (s)W 1 (k,s). (33) 

For Levy walks with 1 < 7 < 2 the center part of the PDF has the form of Levy density 

H3i lag 

W ^~(K^'-'{(K^)' (34) 

where 

1 r+00 f 

with the characteristic function {K 1 is some constant) 

l 7 (q) = exp(-ir 7 |g| 7 ), (36) 

provided that the density was initially localized at x = 0, and exhibits a sharp cutoff 
marked by the ballistic peaks at |x| — vt (17] (clearly seen in figure Ej). 

Taking k = 0, or integrating the PDF equation (133]) in x from zero to infinity, we 
find the average occupation fraction in x > 0, V + (s) = Y + (s)/s. On the other hand, 
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using equation ( 13 Up the average occupation fraction in the material x > in Laplace 
space behaves as V + {s) oc s _1+ ( 1_Q )/ 2 for s — > 0. Comparing two expressions we find 

Y+(s) oc s (1 - Q)/2 , s -> 0, Y+(t) oc r i-(i-«)/2 ) t (37) 

Inverting the Fourier-Laplace transforms of equation ( 133|) and using equations (134"]) . ( 137|) . 
the PDF for x > reads 

P(x, oc ^ dr* ^ / 7 i^-^j , x > 0. (38) 

From equations ( 1381) . ( 1321) it follows that the PDF of particle position for coupled 
sub-super diffusive systems with 1 < 7 < 2 and < a < 1 possesses the scaling form 

jt a/2 P(x,t), forx<0 
Ga ' 7(0 K { t 1 /7 + i/2- Q / 2p(X; t)) for x > Q) ( 39 ) 

and 

f |x|/t a/2 , for x < 

* = Ah * n ^ 4 °) 

[ x/t /7 , for x > 0. 

To reconfirm this result notice that using the scaling of the PDF, the occupation fraction 
in x > is 

V + (t) = jT dx P(x, t) ~ r V7-1/2W2 jf 00 rfx Gq 7 ( JL) „ (41) 

/*oo 

~t- l l 2+a ' 2 ^G Q7 (Ooct- 1/2+a/2 , 

which is what we have found from the FPT analysis (see equation (1301 ). Numerically 
calculated scaled PDF is shown in figure [7] for a = 0.5 and 7 = 1.75. Similarly, for the 
case (II) < a < 7 < 1 

jt a/2 P(x,t), forx<0 
Ga '" i0 K \ t 1+ ^/ 2 P(x, t), for x > 0, (42) 

and 

f |x|/t a/2 , for x < 

\ 11 f n ^3) 

I x/t, tor x > U. 

For the case (III) < 7 < a < 1 the scaling of the PDF for x > is determined by the 
ballistic regime of Levy walks (|x|) oc t [19], 120] 

( t a -^ 2 P(x,t), forx<0 

QO « , p , . f ^ n ( 44 ) 

I i P(x, t), for x > 0, 

and 

' \x\/t a/2 , for x < 

f n (45) 

x/t, tor x > 0. 
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Figure 8. Simulation of the drift for coupled sub-superdiffusive systems with 

q + = q~ (symbols) are favorably compared with our theory equation (|50j) (dashed, 
dotted and dashed-dotted lines). (I) < a < 1, 1 < 7 < 2: a = 0.3, 7 = 1.7; (II) 
< a < 7 < 1: a = 0.35, 7 = 0.8; (III) < 7 < a < 1: a = 0.7, 7 = 0.3. 



3.3. Drift in coupled sub-superdiffusive system 

Using the scaling form of the PDF it is easy to estimate the sign and the time dependence 
of the the mean position of the packet, which is initially at x = 0. As it follows 
from equations (I40j) . (H3l) and (I45j) . Levy walks always spread further than subdiffusive 
trajectories: In all cases for x < 0, \x\ oc W 2 while for x > 0, x oc ^Vt+V 2- "/ 2 f or 
case (I) and x oc t for cases (II) and (III). Therefore, the sign of the mean is always 
positive, that is the drift is directed to the Levy walk independently of q + /q~ . 
Now we calculate the time dependence of the drift 

/oo rO poo 

dx x P(x,t) = / dx x P(x,t) + / dx x P(x,t). (46) 
-00 J — 00 Jo 

For the case (I)0<a<l,l<7<2, using the scaling form of the PDF equations ( |39l 
M), we find 

f° dx x P(x, t) ~ r a/2 f° dx x G an ( ] ~ (47) 
j— 00 J— 00 \t a ' J 

~ e' 2 f di z c an (0 ~ t a/2 , 

J™ dx x P(x, t) ~ t-i/7-i/2+ a /2 dx x Ga ^ „ (48) 

/•oo 

~ t 1 /7-V2+a/2 / Gq (£) _ t l/ 7 -l/2+a/2_ 

j0 
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For 1 < 7 < 2 and < a < 1 the exponent in equation (j48p is always greater than the 
exponent in ( HTj) . 1/7 — 1/2 + a/2 > a/2, and therefor for the case (I) 

(X(t)} OC £ 1 /7-l/2+«/2 > 0) t ^ OO. (49) 

Using equations (I42T) . (t43l) . (t44l) . (|45|) we find the drift for the cases (II), (III). Summarizing 
results we have 

f t^- 1/2+a/2 , (I): 1<7<2, < a < 1, 
(x{t)) oc J t l ^ /2+a/2 , (II): < a < 7 < 1, (50) 

[t, (III): < 7 < a < 1. 

Figure [8] demonstrates a good agreement between the theory and numerical simulations 
for all three cases. Results for the occupation fractions and the drift suggest that for 
the case (I) and (III) (0 < a < 1, 1<7<2 and < a < 7 < 1) the drift is opposite to 
the flow, which for these cases is directed to the subdiffusive part, V~ — > 1 as t — > 00 
(see equation ( 131~|) ). 

Summary 

We investigated systems consisting of two materials with sub or superdiffusive properties 
and a boundary between them. In coupled subdiffusive system particles flow to the 
slower medium, while the direction of the averaged drift is determined by symmetry 
breaking at the boundary, q L 7^ q R in our model. This leads to interesting phenomena 
unique to subdiffusion: (i) under certain conditions all particles are found in one sample 
(e.g. V~ — > 1), but the drift is oppositely directed ((x) > 0), (ii) the drift does not 
depend on properties of the fast medium, namely even if {x(t)) > 0, the anomalous 
diffusion exponent a + and K + in x > do not control (x(t)) (under certain conditions). 
We find similar behavior for the diffusion in quenched trap model and two dimensional 
comb structure, which points out to a broader generality of our results (see [30]). We 
argue that such a behavior is a general feature of subdiffusion in disordered systems. 
For coupled sub-superdiffusive systems we find a net drift, which is always directed to 
the superdiffusive material independently of the asymmetry at the boundary, while the 
direction of the particles flow depends on the relation between sub and superdiffusion 
exponents. These phenomena are explained by the competition of the diffusion processes 
which are slower or faster than normal spreading. 
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Appendix A. Calculation of (n z (t)) 

Here we calculate the drift for the system of two coupled subdiffusive materials. For 
this we need to calculate the average number of returns to the origin, (n z (t)), which is 
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determined in the following way: We define a three state process £(£) = (state 0) if 
the particle is on the origin, = +1 (state +) if the particle is in x > and = —1 
(state — ) if the particle is in x < (see figure iBTl) . In the long time limit the number of 
visits to the origin is independent of Rq since the average waiting times in state + and 
— are infinite. The waiting times in state + are the first passage times, from x = a to 
x — 0, and similarly the waiting times in state — are the first passage time from —a to 
0. These first passage times in the continuum limit were obtained previously [22J and 
they are one sided Levy distributions whose Laplace transform is 



72 



/VK- 



72 



(A.l) 

which implies 4>~(t) oc t _ ( 1+a / 2 ) and similarly for + (t). For a - = a + = 1 we get well 
known distribution of the first passage time of a Brownian motion in half space 
The Laplace transform of the probability to have exactly n z transitions is given by 

z. , s 1 — 4>(s) ~ n Z , \ 



(A.2) 



where <fi(s) = q 4> (s) + q + (fi + (s). Using equation (IA.2j) . the average number of 
transitions to zero is given by 



<«»(»)) 







S i - 018. 



(A.3) 



The small s expansion yields the long t behavior of (n z (t)). Using equations ( lA.lt 



(n z (s)) 



aq^ 



~ < 



VK- 
aq~ 



,-(l+a+/2) 



-(1+0,-/2) 



(q-as a I 2 q + as a+/2S 
S I ; + 



for a > a 
for a~ < a 
for a~ = a 



(A.4) 



Inverting equation (1A.4|) to the time domain and using equation fl5]) we find the drift in 
the long time limit. 



Appendix B. CTRW: Occupation fractions 

Here we calculate the PDF of the occupation fraction in the state — (calculation for the 
state + is similar) for the CTRW. For that aim we use the three state process defined in 
Appendix A. The total time a particle was in the state — after n steps can be written 
as (see figure IBlj) 

n 

t- = J2 T idi + {t-t n )6{*), (B.l) 
i=i 

where 9{ is a random variable with the PDF 

p{0) = q + 6{0) + q-5{0-l), (B.2) 
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Figure Bl. The three state CTRW model constructed from the random walk 
illustrated in figure [TJ 



that is 6i = when the position of the particle is x > or in state + and Oi — 1 
for x < (state — ). In equation (IB. 1[) #(*) denotes the last step and t n = Yli=i r i — 
Yh=i T i@i + Yh=i T i\@i ~ M- The PDF of the occupation time t~ after time t and n steps 
is defined as 

.MO = (6 (r - J2 nOi -{t- t n )0(*)) i(t n < t < t n+1 )\ , (B.3) 



where 



, 1 if condition in parentheses is true, 
I(t n < t < t n+x ) = { (B.4) 
otherwise. 

Now we consider the double Laplace transform of equation ( IB. 31) 

/*oo /*oo 

fs, n (u)= / dte ~ St / dt~ e" ut ~f n 
Jo Jo 



I n,t~ 



'jT tft e- s */(t n < t < t n+1 ) e --Er=i T ->^-«(*^) e (*) ^ . (B.5) 
Averaging over the last step #(*), we get 



I -{s+u)t n _ -{s+u)t n+1 \ 

+q-{- e - u ^i=i Ti6i+u ^i=i n ) . (B.6) 

\ s + u I 

Averaging now over 9i and summing over all jumps we get the PDF of t~ 

oo 

fs(u) = fs,n( U ) = 
n=0 

V 1 -»''" + '» + ^ 1 -»»Vb.7) 



I — (q~ (J)~ (u + s) + q + <fi + (s)) \ u + s 
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where =b are the waiting time PDFs in states + and — . Using the long time limit (or 
small s) ± given by </r(s) oc 1 - as a ~/ 2 /VK- and 0+(s) oc 1 - as a+/2 /\fK+ [22] as 
s — > (see equation ( lA.lj) ). we finally obtain equation (jSJ). 

For q~ = q + and a~ = a + = a equation (jSJ) is reduced to the Lamperti PDF 

[MIES] 



which is a generalization of well-known arcsine law [48J. The method of inversion of 
equation (1B.8|) to time domain is given in 



Appendix C. Remark on the solution of equation (1141) 

We note that the solution of fractional equations (1T4T) and ( 1201) must be used with care. 
While the PDF of particle's position P(x,t) calculated numerically by CTRW model 
perfectly agrees with analytical theory including the jump at the boundary (figure H]) 
and while this solution gives the correct asymptotic behavior of the occupation fraction 
V~ — > 1 (when a~ < a + ) (see equation (fT2]) ). using equations (fill) and (|20|) the 
occupation fraction V~(t) = j_ QO dx x P(x,t) in sample x < within the fractional 
framework is 



r(t)«i-^ / r +v (c.i) 

q-y/W r(l + ^ + ' 



Note that while equation (IC.lj) gives correct leading term, the correction term differs 



from the exact CTRW result ( fl2l) (compare the Gamma functions). Figure [3] illustrates 
the difference of two solutions. Thus, fractional equation works in the long time limit 
and already the first correction to asymptotic solution shows deviation from exact result. 

An alternative to (j2ip boundary condition can be derived by requiring the equality 
of occupation fractions calculated by the CTRW model ( IC.ip with the occupation 
fractions obtained from the fractional equation ( fl2l) 

g q+VWs^^C-is) = q-V~K+s a ~ /2 C + {s), (C.2) 

where Q = 1 + " ~ a . This boundary condition leads to the solution 

C~(s) = %- + - , C + (s) = 2 I _ + . (C.3) 



q- V K+" " ' ^ q+ V K~ 

Analytical probability density function equations ( Tl4]) and ( 1C.3I) is different compared 
to equations ( IT4|) and (T20|) . As the latter was derived from the equality of occupation 
fractions calculated by the CTRW model ( 1C.ll) and the drift obtained with the fractional 
equation (lT2~j) . it does not describe the jump of the PDF at the boundary. This solution 
also gives different results for the moments including the drift. 
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